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We present numerical methods to solve the Generalized Hartree-Fock theory for fermionic sys- 
tems in lattices, both in thermal equilibrium and out of equilibrium. Specifically, we show how 
to determine the covariance matrix corresponding to the Fermionic Gaussian state that optimally 
approximates the quantum state of the fermions. The methods apply to relatively large systems, 
since their complexity only scales quadratically with the number of lattice sites. Moreover, they are 
specially suited to describe inhomogenous systems, as those typically found in recent experiments 
with atoms in optical lattices, at least in the weak interaction regime. As a benchmark, we have 
applied them to the two-dimensional Hubbard model on a 10 x 10 lattice with and without an 
external confinement. 
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I. INTRODUCTION 

Recent experiments with cold atomic systems have 
attracted the interest of several scientific communities. 
Many physical phenomena, like Bose-Einstein condensa- 
tion, superfluidity, or the formation of Cooper pairs, have 
been observed in seminal experiments and carefully an- 
alyzed and characterized [H [2]- In the last years, the 
possibility of loading the atoms in optical lattices offers 
a new avenue to the observation of many other intriguing 
phenomena [3 . 

There exists several theoretical tools in order to de- 
scribe most of those phenomena. For the case of Bosons, 
the Gross-Pitaevskii equation provides us with the basic 
tool to analyze many static and dynamic properties for 
cold atomic gases [4]. Such an equation is valid when- 
ever the occupation number of a particular mode is large, 
and it is equivalent to a mean-field theory. Typically, it 
can be solved numerically, offering very accurate predic- 
tions and explanations of all the observed phenomena in 
the weak interacting regime. Furthermore, it can also be 
used to describe lattice systems, at least as long as strong 
correlations are absent. In the case of Fermions, the basic 
tools at hand are Hartree-Fock and BCS theory [5 . The 
first one assumes that the state of the particles can be 
approximated in terms of a Slater determinant, and thus 
it also correspond to a mean-field theory. The latter as- 
sumes a specific form of the wavefunction, whereby pair- 
ing between Fermions is allowed. For both approaches, a 
time-dependent version has been constructed [6HTQ]. 

Hartree-Fock and BCS theories can be unified in terms 
of a more general framework, the generalized Hartree- 
Fock Theory (gHFT) [11 . The main idea is that both, 
the Slater-Determinant- and BCS-states are members of 
a larger family of states, the so-called Fermionic Gaus- 
sian states (FGS). Therefore, one may try to find the 
state within this family which best approaches the real 
state of the system under study. FGS are those states 
whose density operator can be expressed as an exponen- 
tial of a quadratic function of canonical creation and an- 
nihilation operators. This property immediately implies 
that they fulfill Wick's theorem: that is, all correlation 



functions can be expressed in terms of the second mo- 
ments of all creation and annihilation operators. Such 
moments are gathered in a matrix, the so-called covari- 
ance matrix (CM) which completely characterizes the 
state [iT. Thus, FGS can be very efficiently described, 
since the number of parameters only scales quadratically 
with the number of available fermionic modes and the ex- 
pectation values of observables can be easily computed. 
For the case of lattices, the number of modes is automat- 
ically finite, and thus it should be an ideal playground 
for such a gHFT. In fact. Bach, Lieb and Solovej ^ 
were the first who applied such a theory to a Hubbard 
model in 2D, and where able to solve the corresponding 
equations exactly for the homogeneous case, both for the 
ground state and in thermal equilibrium. 

gHFT can be easily applied to any lattice system, ho- 
mogeneous or not. The resulting equations are, however, 
difficult to solve in general. In this paper we derive such 
equations for the CM, and introduce different methods 
to solve them. First of all, we derive the equation for the 
real time evolution of the CM. Second, we consider the 
state within the FGS which minimizes the total energy, 
and thus approximates the ground state. For that, we 
derive the evolution equation (in imaginary time) for the 
CM in the presence of a lattice Hamiltonian. We show 
that under such an equation, the CM converges to the one 
which minimizes^ the energy within the family of FGS. 
Thus, the evolution equation provides us with a practi- 
cal way of determining the approximation to the ground 
state. Finally, we derive the equation for the CM that 
minimizes the free energy, and apply a simple fixed-point 
iteration method to solve it. 

Then we apply the methods to the spin- 1/2 Hubbard 
model in a 10 X 10 lattice. First, we benchmark them 
with exact solutions of the gHFT (for the homogeneous 
case) [11], finding an extraordinary agreement both for 
attractive and repulsive interactions. Second, in order to 



^ Strictly speaking, the CM converges to the solution for which 
the energy is extremal. 



test the validity of gHFT itself, we compare the results 
for positive interactions, where strongly correlated effects 
are expected to be most pronounced, with the Monte 
Carlo results of Refs. [T3", IT. The results are obviously 
less precise than those obtained with such methods, other 
Monte Carlo approaches [15 , or the recently developed 
algorithms based on PEPS and other tensor networks 
states [16-21 . Although we obtain a good qualitative 
agreement with [13l [14] in many regimes, the FGS are, 
as expected, not able to capture all possible fermionic 
phases in the strong-correlation regimes, as they consti- 
tute a subclass of all possible fermionic states. For ex- 
ample, the finite temperature Mott phase cannot be ap- 
proximated well as a Gaussian state, and is thus absent 
in the phase diagram. 

At this point we would like to warn the reader that in 
this work we will mainly concentrate on the possibilities 
offered by the numerical approaches, and not focus on the 
physics of the Hubbard model. To be precise, our goal 
is a) to formulate the standard gHFT for fermionic sys- 
tems with a two-body interaction in the lattice using the 
language of covariance matrices, b) to derive numerical 
methods to determine the ground and thermal states as 
well as the time evolution of these systems within gHFT 
and c) to benchmark our techniques by applying them to 
the two-dimensional Hubbard model. 

This paper is organized as follows: In Sec. [H] we give a 
precise formulation of the problem we aim to solve, and 
introduce the concepts necessary for the understanding 
of FGS and gHFT. This Section includes further a dis- 
cussion why FGS are supposed to capture well the phys- 
ical properties of some of the most relevant fermionic 
models. Next, we derive the methods to approximate 
the real time evolution, the ground and thermal states 
of fermionic systems within gHFT in Sees. [Ill A - IIIC 



Then, we benchmark our methods by applying them to 
the two-dimensional translationally invariant Hubbard 
model, both in the attractive and the repulsive regime, 
in Sec. |IV| In addition, motivated by the current possi- 
bility of implementing the Hubbard model in more than 
one dimension with the help of optical lattices, we also 
investigate the 2d Hubbard model with an external con- 
finement in Sec. [V| and close with a study of dynamical 
processes, both for the translationally invariant and the 
trapped system, in Sec. |VI[ 



II. STATEMENT OF THE PROBLEM 

In this Section we give a statement of the problem, 
summarize the properties of fermionic Gaussian states 
that are necessary for the understanding of our work (for 
more details see e.g. [12 ) and argue why these states 
capture the properties of many relevant physical models. 

In our work we aim at understanding the physical prop- 
erties of fermionic systems in a lattice, like for instance 
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where the ak denote fermionic mode operators obeying 
canonical anti-commutation relations (CAR), {a^^a]} = 
5ki and tkuUkimn ^ C. These models describe fermionic 
atoms in optical lattices [22 , and one of the most promi- 
nent ones, the Hubbard model ^, is expected to de- 
scribe high-Tc superconductivity [24]. We approach the 
problem in the Major ana picture and define Ck = a\^ak^ 

Ck-\-M = (— 0(^1 ~ ^k)i where M is the number of modes 
and k = 1,...,M. These operators satisfy the CAR 
{c/c,q} = 2Ski. The Hamiltonian H in this picture is 
given by 



H(T,U) = i^^TklCkCi + ^ UklmnCkClCr, 



(2) 
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where TkuUkimn ^ K.- The CAR allow us to antisym- 
metrize T and U such that T^ = — T while U is an- 
tisymmetric under the exchange of any of two adjacent 
indices. 

Our goal is to derive ground and thermal states as 
well as the time-evolution of systems characterized by 
the Hamiltonian given in Eq. ([l]) within the family of 
FGS. These are states that can be represented as an ex- 
ponential of a quadratic form in the Majorana operators, 
and thus they are fully characterized by their second mo- 
ments collected in the real and anti-symmetric covariance 
matrix (CM) F^^ = (f [c/c,q]) from which all higher cor- 
relations can be obtained via Wick's theorem (see e.g. 

iPtr[pc,-,...c,,J=Pf(r|,,...,.,J, (3) 

where 1 < ji < . . . < J2p < 2M and rj^^...^J2p is the 
corresponding 2p x 2p submatrix of F. Pf(rj^^...^J2p)^ = 
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J2p. 



is called the Pfaffian^. F is the CM of a 



physical state iff zF — 1 < 0, while pure states have to 
fulfill F^ = —1. Every pure FGS is the ground state of a 
quadratic Hamiltonian 
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with real and antisymmetric Hamiltonian matrix h. All 
Gaussian states remain Gaussian under the time evolu- 
tion governed by a quadratic Hamiltonian, and the CM 
transforms according to T{t) = O(t)F(0)O(t)-^, where 



0{t) 



^4ht 



is an orthogonal transformation. 



Let us now explain why FGS provide a very pow- 
erful technique for the description of fermionic many- 
body systems. First, note that all pure states within 
this family can be brought into a standard form |^) = 
Yli.{uk + v/c«lttL/e)|0), where \uk\'^ + \vk\'^ = 1 V/c via a 
change of basis [25 . Thus, the pure Gaussian states in- 
clude the BCS-states introduced in the theory of super- 
conductivity, as well as all Hartree Fock states. Second, 



The sign of the Pfaffian is determined by the condition that the 
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appears with positive sign. 
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FIG. 1: a) The time evolution under the interaction Hamilto- 
nian M within the set of Gaussian states (purple) can be re- 
formulated as the time-evolution under a quadratic but state- 
dependent Hamiltonian Heff = '^^ijhij{T)ciCj, where h is 
given in Eq. (Is]), that respects the Gaussian structure, b) 
The minimization of the energy within the set of Gaussian 
states can be mapped to an imaginary time evolution within 
the set of Gaussian states with the effective, state- dependent 
Hamiltonian given by Heff — i X^- • hij{T)ciCj^ with h defined 
in Eq. (|8|. c) The minimization of the free energy leads to 
the implicit equation F = ztanh [2i;^/i(F)] (c.f. (26)) for the 
covariance matrix, that can be solved, e.g., by a fixed-point 
iteration. 



FGS also include all mixed states for which Wick's theo- 
rem ([3| applies, e.g., thermal states of quadratic Hamil- 
tonians. Thus, all together, we see that FGS include a 
variety of states which form the basis of several many- 
body phenomena. 



III. INTERACTING FERMIONS IN 
GENERALIZED HARTREE FOCK THEORY 

In the following Section we derive the theoretical 
framework that will allow us later on to simulate 
fermionic many-body problems, like the Hubbard model, 
for big system sizes. This is achieved by approximating 
the state of the fermions by a FGS, something which is 
implemented by means of Wick's theorem. 



A. Real Time Evolution 

We start the section deriving an evolution equation 
that allows us to describe the time-evolution of a system 
at zero and finite temperature evolving under the inter- 
action Hamiltonian (fTj) within generalized Hartree-Fock 



theory. The key idea is to use the covariance matrix for 
the description of the dynamical evolution. In the Heisen- 
berg picture, where the time evolution of the Majorana 
operators is given by Ck{t) = e'^^i^^u)t^^^-iH{T,u)t ^ ^^^ 

covariance matrix evolves according to 



d 
dt 



T^^{t) = {[ca,{t)c^{t),H{T,U)]). 



(5) 



As H{T^ U) involves terms quartic in the Majorana oper- 
ators this evolution clearly takes us out of the Gaussian 
setting. We truncate this transformation imposing that 
Wick's theorem holds, i.e. we write 

{ciCjCkCi) = -{TijTki - TikTji + TiiTjk). (6) 

With the help of the commutation relation 

[CcC/3, CiCj] = 2{CiCaSjf3 " CiCfsSja " CfsCjSia + CaCjSi^), 



we can calculate the contributions for the quadratic 
term Hq = Y.ki ' 



Lq — '}2m TkiCkCi and the interaction term Hj 



{[caC0,Hq]) = 4[r,r]„^, 

This imphes the fohowing time evolution of the covari- 
ance matrix: 
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where we have defined tTB[UT]ij = "^i^iUijki^ik- This 
equation can be formally integrated: 



Tit) = o{t)r{o)o{tf, 



0{t) 



= Texp ^4 / h{T{t'))dt'\ , 



(9) 



where T denotes the time ordering operator. Note that 
due to the anti- symmetry of h{r{t')) Eq. ([9| guaran- 
tees that the matrix 0{t) is an orthogonal transforma- 
tion. Hence, when starting with a valid covariance matrix 
r(0), this approximation scheme ensures that we remain 
within the set of Gaussian states. Further, it follows from 



what we have said in Sec. HI that the dynamical evolu- 
tion under the interaction Hamiltonian H(T^ U) can be 
understood as the time evolution under a quadratic hut 
state- dependent Hamiltonian 



Hq{V) 



i^h{V{t))kiCkCi. 
ki 
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Further, this approximation scheme does not only ensure 
that we always remain within the set of Gaussian states, 
but a short calculation shows that it also conserves energy 
and, for a number conserving Hamiltonian, the particle 
number. We start with the energy, which is given by 



E{t) = tT[Hp] = -tr[(T + 3tTB[UT])r], and hence, with 
the help of Eq. ([7|), it fohows that 



dE 



^W = E ar^r,, = 4tT[h{T)[h{r),r]] = o. (ii) 
ki ^^ 

To prove the conservation of the mean particle number 
A^, as long as [H^ TV] = 0, we write the particle number 
operator N = Y.k ^l^k = T" + i ^k,i ^kiCk^u where v is 
real and anti-symmetric, and obtain N{t) = ^ — |tr[z/r]. 
Thus (c.f. Eq. ^) 

j^N{t) = -itr[^f] = -i^^,,tr[p[QC,,i7]] 

kl 

= -itT[p[N,H]] =0. 

Hence, we have shown that applying Wick's theorem 
to the evolution equation of a system governed by a two- 
body interaction Hamiltonian leads to a consistent dy- 
namical equation for the CM. The truncation of the evo- 
lution is equivalent to an evolution under an effective 
state-dependent quadratic Hamiltonian within the set of 
Gaussian states, as depicted in Fig. [l]a). We remark 
that this approach works for pure as well as mixed FGS. 
Note also that the remaining equation resembles in some 
way the Gross-Pitaevskii equation for Bosons since that 
one can also be interpreted as giving the evolution in 
terms of a Hamiltonian that depends on the state; how- 
ever, the main difference is that in the case of Bosons 
such an equation is obtained for the first moments of 
the field operators. Note that in the case of Bosons one 
could also perform the same approximation for the first 
and second moments, obtaining a set of coupled equa- 
tions which would include the condensate part as well as 
thermal particles. 



B. Ground States 

In principle, the ground state in gHFT can be found 
via a direct minimization 
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(12) 



For the translationally invariant Hubbard model this 
problem could be reduced to an optimization over two 
parameters only [11]. Generically, this constrained 
quadratic minimization problem is a daunting task. We 
attack this problem from a different perspective. A well- 
known approach for the determination of the ground 
state l^o) of a Hamiltonian H is imaginary time evo- 
lution. Due to the exponential growth of the state space 
with the number of modes this approach can be applied 
to small systems only. The idea is to apply the Gaussian 



approximation (in the form of Wick's theorem) to derive 
an evolution equation for the CM. In this way we get a 
simple quadratic scaling in the computational time with 
respect to the number of modes. 

Naively, one would think of obtaining the evolution 
equation in imaginary time by replacing t \-^ it in the 
real-time evolution equation of the CM, Eq. ^. How- 
ever, as it turns out, this is not the right road to take. 
Though, as we will show below, the following approach 
will lead us to the desired ground state: We start with 
an arbitrary pure Gaussian state p(0), and evolve it ac- 
cording to 
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where the normalization of the density operator for any 
time, i.e. the denominator tr[e~^^^p(0)], is crucial for 
obtaining the correct equations. Since we apply the non- 
quadratic operator H to the Gaussian state p, this proce- 
dure clearly takes us out of the setting of Gaussian states. 
One way around would consists of discretizing the evolu- 
tion and finding for small time steps At the best Gaussian 
approximation at each step. However, due to the trunca- 
tion it is not clear that this procedure will converge. And 
even if we find a steady state it is not clear that this will 
be the ground state of the system. Another possible ap- 
proach is to use the quadratic but state-dependent effec- 
tive Hamiltonian HgiT) that is derived from H(T^ U) for 
the imaginary time evolution. As Hq{T) is a quadratic 
operator we will stay in the space of Gaussian states. 
But as the Hamiltonian is state-dependent the outcome 
of this procedure is not clear. However, as we will show 
below, both approaches are equivalent and lead indeed 
to the desired solution. To be precise, we show that the 
following is equivalent: 



1. Direct minimization of the energy Eq. (12) in gen- 
eralized Hartree-Fock theory. 

2. Imaginary time evolution of the state p with the 
full Hamiltonian H{T^ U) for small time steps At 
followed by an approximation of p{t + At) by a 
Gaussian state. 

3. Imaginary time evolution of p with the quadratic 
but state-dependent Hamiltonian Hq{T). 



1 . Minimization of the energy 

In order to obtain the generalized Hartree-Fock ground 
state we have to solve the minimization problem 



min E{p) = min ii:[Hp\ 

p Gaussian p Gaussian 



]^^, { E^^i^i,- - 3 ^ UijkiTijTu \ . (14) 

i,j,k,l 



ir<i 



It has been proven in Ref. [TT] that the HF ground state 
is always pure, i.e. F^ = — 1. Using Lagrange multipUers, 
we arrive at the fohowing necessary conditions for a local 
minimum (see Appendix H)^: 
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These two equations are non-linear matrix equations and 
thus hard to solve, both analytically and numerically, for 
large systems. But we show next that these equations 
appear as the steady-state conditions of imaginary time 
evolution. 



2. Imaginary time evolution 



From Eq. (13) we see that the evolution of the density 



operator p under any Hamiltonian H in imaginary time 
is given by 

m = -{H,p{t)} + 2p{t)tr[Hp{t)], (17) 

SO that the covariance matrix evolves according to 

d 



dt 



Tkiit) = -itT[{H, CkCi}p{t)] + 2r kitT[Hp{t)]. (18) 



We show in Appendix [B] that both approaches for imag- 
inary time evolution, number 2 and 3, lead to the same 
evolution equation of the covariance matrix: 



d 
dt 



F = -4 (F/i(F)F + /i(F)) 



(19) 



First, note that Eq. ([19| ensures that a pure state re- 
mains pure under the evolution, since 



d 
dt 



F^ = Ff + f F = -4(F^/iF + F/i + F/iF^ + hT) = 0, 



as all pure states fulfill F^ = —1. Next, we show that the 
energy always decreases under the evolution: 



d BW - 

^W = E^T^r,,=^MF),,f,, 

ki ^^^^ ki 

= ^ir[h{VhV^h)] =2ir[{[hJ]f] < 0. (20) 



dt 



Here, we have used that the matrix [/i, F] is antisymmet- 
ric, implying that ([/i, F])^ is negative-definite. Further, 
it follows from Eq. ([20| that the steady-state is reached 
iff [/?',F] = 0. Thus, we are ensured to reach the gHF 



^ Note that Eqs. \lb\ and ( |16| describe all minima, local and 
global ones. However, for the particular problems studied in 
the next sections, our numerical investigations show that our 
algorithms lead in general to the global minimum, i.e. the ground 
state. 



ground state via an evolution with Eq. ( 19 ) with initial 
condition F(0)^ = -1. 

Finally, if we start in a pure input state, i.e. F(0)^ = 
— 1, Eq. ([l9| can be formally integrated (see also 

Fig.[l}D): 



0(t)T{^)0{tf 
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(21) 
(22) 



The evolution of F with the orthogonal 0{t) has turned 
out to be especially fast and stable in numerical imple- 
mentations compared to the direct integration of Eq. 
(19) via e.g. a Runge-Kutta solver, and we have used 
it for our numerical applications. 



C. Thermal states 

As we have seen in the last Section, the generalized 
Hartree-Fock ground state can be obtained via an imag- 
inary time evolution. To learn about the finite tempera- 
ture properties in gHFT we have to consider the approx- 
imation to the Gibbs state p ~ e~^^, where f3 = -^^ is 
the inverse temperature. We recall that the Gibbs state 
minimizes the free energy, F{p) = E{p) — j3~^S. Thus, 
we have to solve 



min F{p) 

p Gaussian 



min {E{p) 

p Gaussian 



r^s{p)) 



As in the case of the minimization of the energy, we re- 
formulate the problem as an optimization problem using 
Lagrange multipliers. As we show in Appendix [C) this 
approach leads to the following necessary conditions for 
a minimum of the free energy: 
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(23) 
(24) 



(25) 



Now, assume for the moment that det(l+F^) ^ 0. Then, 
according to Eq. (24), /i^? = is the only possible solu- 
tion. Hence, the CM has to be of the form 



itanh [2if3h{T)] . 



(26) 



Now, a CM of this form always fulfills Eqs. ([23|) and ( |24[ ), 
and includes, in the limit /3 ^ oo, the zero-temperature 
case, so that Eq. (26) is indeed a solution of Eqs. (23) 
and (l24l). 



Note that from the perspective of numerical implemen- 
tation, we can solve Eq. (26) e.g. via a fixed-point iter- 
ation (see Fig. lip). 
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FIG. 2: Ground state ener gy o f the translationally invariant 
Hubbard model (c.f. Eq. ( [27| ) for an interaction strength 
u G [—10, 10] at half-filhng, i.e. /j = 0. The blue line shows 
the exact gHF ground state energy, ^exact, obtained from 
Ref. [11], while the triangles correspond to the energy ob- 
tained via imaginary time evolution, Eimag- 



FIG. 3: Energy as a function of the chemical potential for the 
attractive Hubbard model at various values of the interaction 
strength u. The dotted lines depict the exact gFH solution, 
while the triangles and diamonds correspond to the solutions 
obtained via our approach. We find excellent agreement of 
the two approaches 



IV. BENCHMARK: THE TWO-DIMENSIONAL 

TRANSLATIONALLY INVARIANT HUBBARD 

MODEL 



In the following, all calculations are performed for a 
10 X 10 lattice. 



In the following Section we apply our method to 
the two-dimensional translationally invariant Hubbard 
model with periodic boundary conditions: 

where x and y are points on a two-dimensional lattice, 
(i,j) denote nearest-neighbors and a =t,i denotes the 
spin degree of freedom. We consider only the case where 
we have the same number of spin-up and spin-down par- 
ticles, so that we can use the same chemical potential for 
the two species. For u < {u > 0) the second term in 
H is an attractive (repulsive) on-site interaction between 
particles of opposite spin. In the following, we set t = 1 
as the energy scale of the system. Note that the case of 
half- filling is characterized by /i = 0. 

The goal of this Section is two- fold: In the first part. 
Sec. IV A[ we benchmark our numerical method by con- 
sidering physical quantities for which an exact solution 
within gHFT is known. In the case of the translation- 
ally invariant Hubbard model it could be proven in Ref. 
[11] that the (free) energy for the (thermal) ground state 
can be found via a two-parameter optimization, and we 
will compare our numerical results with the numbers ob- 
tained from the optimization. After the demonstration 
of the power of our approach, we continue in the second 
part, Sec. |IVBl with a comparison of gHFT itself to more 
powerful and and sophisticated methods, like Quantum 
Monte Carlo (QMC). 



A. Comparison with the exact gHF solution 

For the translationally invariant Hubbard model 
Eq. ( [27| ) the exact results for the energy of ground and 
thermal state within gHFT where presented in [11]. In 
Fig. [2] we compare the ground state energy for half- 
filling obtained via imaginary time evolution accord- 
ing to Eq. (21) (triangles) with the exact gHF solu- 
tion (blue line) and find excellent agreement. The rel- 
ative error in the energy is given by ((i£^)rei = |£^imag — 

^exactl/l^exactl < 2 • 10~^. 

In Fig. [3] we show results away from half-filling. In 
this case, a closed solution for the energy could only be 
provided for the attractive Hubbard model in Ref. [TT] . 
We have compared the exact gHF solution (dotted lines) 
with the results obtained via our approach (triangles and 
diamonds) for u = —2, —6, —10, proving that our numer- 
ical method also works well for a doped system. 

Next, we consider the case of finite temperature, where 
the gHF Gibbs state is given by the implicit equa- 
tion (26). Starting from the ground state we change 
the temperature in steps A/3 = 0.01 and compute the 
new thermal state via a fixed-point iteration. In Fig. [4| 
we present a comparison of our solutions (triangles and 
diamonds) for half-fiiling with the exact results (lines) 
for various values of ix in a temperature range of /3 G 
[0.4, 1.6], and find excellent agreement. 

One of the most interesting questions concerning the 
Hubbard model is the appearance of superconductivity. 
As we have already pointed out in Sec. [H] FGS describe 
among others superfluid and normal states, since every 
pure state of this family can be brought into the standard 
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responding transition temperature has been derived. In 
Fig. |5] we have depicted the pairing as a function of the 
inverse temperature for various values of the interaction 
u^ and find a breakdown at a finite value Pc that depends 
on the interaction strength. The values of /3c that we ob- 
tain via our numerical method agree well with the results 
presented in Ref. [Tl], and that are depicted as circles in 
Fig. [5J Further, we give a list of the critical exponent 
7 defined via P = a{Tc — T)^ for various values of u in 
table [D 
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FIG. 4: Free energy of the translationally invariant Hubbard 
model (c.f. Eq. (27)) in the temperature range /3 G [0.4, 1.6] 
for iz = -1, . . . , -10 at half-filling. The lines depict the gHFT 
solution given in Ref. 11 , while the triangles and diamonds 
show the results obtained via the fixed-point iteration. 



TABLE I: Critical exponent 7 for the phase transition from a 
paired to an unpaired phase for various values of u. We have 
fitted P around the critical temperature Tc to a curve of the 
form P = a{Tc-Ty for lattice sizes 8 x 8, 10 x 10 and 12 x 12, 
finding the same values of 7 for the given precision. 




FIG. 5: Temperature-dependence of the pairing (see Eq. ( p8| ) 
for various values of the interaction u at half-filling. We ob- 
serve the occurrence of a phase transition from a paired phase 
to an unpaired normal state. The circles indicate the value of 
the transition temperature derived for gHFT in Ref. [IJJ. 



(BCS) form |^) = Ylki'^k + ^/c<^IttL/e)|0) via a change of 
basis. In this basis, the gap parameter A = ^^ uj^Vk 
is an order parameter for the superfluid phase. To con- 
sider superfluidity independent of the basis, it is more 
appropriate to consider the normalized pairing measure 
derived in [26]: 



M 



Ek4«') 



(28) 



kl 



Then, a positive value of P indicates that we are in a 
paired phase, while a normal state is characterized by 
P = 0. 

It has been shown in Ref. \TV that within gHFT the at- 
tractive Hubbard model exhibits a phase transition from 
a paired phase to an unpaired normal phase at a critical 
temperature Pc for any flnite lattice size, and the cor- 



Finally, we give in Fig. [6] a phase diagram of the pair- 
ing for the ground state as a function of the chemical 
potential and the interaction strength. (Note that the 
value of the pairing is not unique in the case of half- 
fllling due to the Gauge symmetries of the Hamiltonian 
(see Ref. [H]). We have restricted to represent the phase 
diagram for negative values of u, since we flnd that P = 
for u positive. 

In summary, we have demonstrated that our numer- 
ical methods are capable of obtaining the correct gHF 
ground and Gibbs state for the translationally invariant 
Hubbard model, both in the attractive as well as in the re- 
pulsive regime. Having demonstrated the validity of our 
approach, the next obvious question is now how gHFT 
itself compares to more powerful and sophisticated meth- 
ods, like Quantum Monte Carlo (QMC). 



B. Comparison with Quantum Monte Carlo 

In this Section we compare the results obtained via our 
algorithm to numerical data from QMC simulations. To 
this end, we will concentrate on the repulsive Hubbard 
model at half-fllling. The reason for this choice is two- 
fold: First, as it has been shown in [271 EH], the physics of 
the attractive Hubbard model is well-described by BCS- 
theory, which is included in gHFT. This motivates to con- 
sider primarily the repulsive Hubbard model, and com- 
pare our results with results obtained from QMC. While 
QMC suffers from the notorious sign problem when ap- 
plied to fermionic systems, it could be shown that in 
the case of a half-fllled lattice the sign problem does not 
occur [29], and thus QMC becomes exact. As a con- 
sequence, the two-dimensional repulsive Hubbard model 
at half-fllling has undergone an intense investigation in 
recent years. A partial list of these results is given in 
j3QH4Qj and references therein. In this Section, we com- 
pare our results obtained in gHFT with the recent QMC 
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FIG. 6: Pairing (see Eq. ( 28 )) of the ground state as a function 



of the interaction strength and the chemical potential in the 
translationally invariant case. As a guide for the eye, we have 
added lines for the particle numbers 50, 100 (half filling) and 
150. 



results obtained in Ref. [13 and Ref. [TT for half- filled 
and doped lattices. 



1. Ground state properties for half-filling 

We start with a comparison of our approach with the 
results presented in Ref. [13] , where a half- filled system is 
considered. Since in the latter work the numerical results 
are presented for various lattice sizes, while our results 
are always given for a 10 x 10 lattice, we have to restrict 
to a qualitative comparison in certain cases. As the first 
physical quantity, we consider the momentum distribu- 
tion n(k) that has been depicted for u = 2,3,4,5,6,8 
and very small, but non-zero temperature in Ref. [T3]. 
In Fig. |8]we show our results for the momentum distri- 
bution n{k) for an interaction u = 1, . . . , 10 at half- filling 
in the entire Brioullin zone. Our results show an excel- 
lent qualitative agreement with Fig. la of Ref. [13 . As in 
Ref [13 we find that with increasing u^ the momentum 
distribution fulfills n(0, 0) ^ 1 and n(7r, 0) = 1/2 for all 
values of the interaction strength, and n{k) grows with 
increasing u for k G [(7r,7r), (tt, 0)]. Further, as in [13] 
we observe a broadening of the Fermi surface with grow- 
ing interaction strength, i.e. the momentum distribution 
gets smeared out. 

In the following, we consider two-particle properties of 
the system. To this end, we study as in Ref. [13] the 
following magnetic properties of the system: 

• The equal-time spin-spin correlation function 

• The local moment (7(0,0) = ((n^^ — n^^)^). 
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FIG. 7: Comparison of the ground state energy obtained 
from gHFT with the QMC results of Refs. [30] and [31]. 
We have depicted the relative error e(u) — \Eqmc{u) — 
EgHFT{u)\/\EQMc{u)\ for lattice sizes 8x8, 10 x 10 and 
12 X 12 and interaction strengths it = 2, 4, 8, 10. 
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FIG. 8: Momentum distribution n(k) of the ground state of 
the repulsive Hubbard model with an interaction ranging from 
u — 1, . . . , 10 at half- filling. While the Fermi surface is sharp 
for a weak coupling, we find a broadening with increasing u. 
(Compare with the QMC result. Fig. la of Ref. 13 .) 



• The magnetic 

E^e^^-^C(f). 



correlation function S{k) 



In Fig. [9] we have depicted the local magnetic moment 
of the ground state, C(0,0) = ((n^^ — n^^)^), for an 
interaction strength u = 0, . . . , 10 at half-filling. The 
purple squares depict the results obtained via our algo- 
rithm, while the blue diamonds represent the QMC re- 
sults of [13]. Singly occupied sites fulfill C(0,0) = 1, 
while for empty or doubly occupied sites C(0, 0) = 0. In 
the non-interacting regime, all four states (full, empty, 
single occupied in the two spin-states) are equally prob- 
able, so that C(0,0) = \. We find that the on-site re- 
pulsion suppresses doubly occupied, and for a half- filled 
lattice also the empty configuration. With increasing u 
the singly occupied configuration is enhanced, leading to 
a well- formed moment at each site. 

Next, we study the equal-time spin-spin correlation 
function C{y) = ((n(^+^)^ - n^£+^)i){n^t - ^xi))- This 
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FIG. 9: Local moment C(0, 0) = {{n^^ - n^^Y) of the 
ground state of the repulsive Hubbard model for an inter- 
action u — 0, . . . , 10 at half-filling. Starting from the non- 
interacting limit, where C(0, 0) = |, the system approaches a 
value of C(0, 0) = 1 with increasing interaction, correspond- 
ing to a phase with no double-occupancy. While the purple 
lines and squares depict the values of the magnetic moment 
obtained with our algorithm, the blue triangles represent the 
QMC results of Ref. [13], Fig. 4. 
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FIG. 10: Equal-time spin-spin correlation function C{y) — 
((n(^+^)^ - n(^+^)^)(n:^^ - n^+i)) of the ground state for an 
interaction u — 1, ... 10 and half- filling. Antiferromagnetic 
correlations increase with the interaction strength, but extend 
even for small interactions over the entire lattice. Compare 
with Ref. 13 Fig. 7, where C{y) has been obtained for a 
24 X 24 lattice. 



FIG. 11: Magnetic correlation function S{k) = Y.T^'^'^C{1) 
for the ground state for u — 1, . . . , 10 at half- filling. The sharp 
peak at (tt, tt) emphasizes the antiferromagnetic correlations 
inherent in a half-filled lattice. (Compare with Fig. 8 of 
Ref. pjj) 
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FIG. 12: Local density rix = rix^ + rtxi as a function of the 
chemical potential for the ground state at n = 6. (Compare 
with Fig. 2a of Ref. [14) 



S{k) is depicted for u = 2). Further, we see that the 
magnetic correlations depend strongly on the interaction 
strength. 



quantity measures the extend to which two spins on site 
X and x^y are aligned. Note that the definition of C{y) 
is rotationally invariant. Our results are presented in Fig. 
[Tol The correlations extend, already for small values of ii, 
over the entire lattice, as it is predicted for a half- filled 
lattice. As one would expect intuitively, the antiferro- 
magnetic correlations are enhanced with growing inter- 
action strength. Our results show the same qualitative 
behavior as in [13], where C{yi) has been calculated for 
a 24 X 24 square lattices. 

Finally, we consider the magnetic correlation function 



S(k) = Xlf^^ ^(0 i^ ^ig- 11 ^ sharp peak occurs 



at momentum (7r,7r), emphasizmg the antiferromagnetic 
correlations. (Compare with Fig. 8 of Ref. [13 , where 



2. Ground state properties of a doped system 

Next, we consider the case of a doped system, i.e. 
/i 7^ 0. This instance has been studied e.g. in Ref. [14]. 
There, the local density Ux was calculated as a function of 
the chemical potential. The results obtained via our ap- 
proach (see Fig. p!2|) show a qualitative agreement with 
Fig. 2a of Ref. |14j, where a density of n{/i) = 1 is 
found for /i ^ — 1 , . . . , 1 . We find that for /i = we have 
half-filling, as expected, while an increase in the chemical 
potential results finally in doubly occupied sites. 

We include here a remark on the limitations of gHFT. 
It is believed that the doped Hubbard model with re- 
pulsive interaction is a highly correlated state, exhibit- 
ing pairing for some value of the doping. However, for 
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FIG. 13: Antiferromagnetic correlation A = {nx^nx-\-di.) in for 
the repulsive Hubbard model at half-filling as a function of 
the inverse temperature /3 G [0, 1] for an interaction strength 
1^ = 10 and lattice spacing d = 1, . . . , 10. We find that the 
antiferromagnetic order is destroyed with increasing temper- 
ature due to thermal fluctuations. 



u G [0, 10] and ja G [—5, 5] we only find an unpaired phase. 
This is in agreement with Thm. 2.11 of Ref. [11 , where 
it has been proven that for any positive definite potential 
U the gHF ground state is unpaired, independent of the 
doping. 



3. Thermal state properties for half-filling 

In this Subsection we consider some finite temperature 
properties of the repulsive Hubbard model, and concen- 
trate again on the case of half- filling. We have seen in 
the last Section that gHFT provides us with a powerful 
tool to study the magnetic properties of fermions with re- 
pulsive interactions in a lattice. Especially, we have seen 
that the ground state of the system shows antiferromag- 
netic correlations, already for small values of the inter- 
action strength. These are expected to decrease with in- 
creasing temperature, since the thermal fluctuations de- 
stroy the antiferromagnetic order. We find indeed this 
behavior, which is depicted in Fig. [13J There, we have 
plotted the order parameter 



A{d) = (nx^rix^di) 



(29) 



as a function of the of the inverse temperature for u = 10. 
We have also derived the critical exponent a given by 
A{T) - (T - Tc)" for lattice sizes 8 x 8 and 10 x 10, 
leading to a ?^ 1.3 — 1.4. 

The antiferromagnetic phase is further characterized 
by the fact that the local fluctuations in the particle num- 
ber vanish, i.e. we are in a Mott phase. This phase is 
predicted to last for even higher temperatures than the 
antiferromagnetic phase ("Mott-gap"), since additional 
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FIG. 14: Value of the antiferromagnetic order parameter 
Ai,2 = ^(1) (c.f. Eq. (29)) of nearest-neighbors (solid pur- 
ple line) and order parameter for the Mott phase, Om — 
{(tlxY) — {rix)'^ (dash-dotted blue line), for the attractive Hub- 
bard model with interaction strength u = 25, 50, 75, 100 as a 
function of the inverse temperature, (3 G [0,0.4]. We find an 
absence of the Mott-gap. 



energy is needed to delocalize the fermions [41 . We in- 
vestigate this behavior in Fig. [14] There, the solid pur- 
ple line shows the value of ^41^2 = ^(1) (c-f- Eq. (29)), 
the antiferromagnetic correlations of nearest-neighbors, 
while the blue dotted line is the order parameter for the 
Mott phase, Om = {{"^xY) — (^x)^- We see that indepen- 
dent of the interaction strength the Mott phase as well 
as the antiferromagnetic phase break down at the same 
temperature, i.e. we do not find a Mott-gap. This obser- 
vation can be explained by the fact that the FGS cannot 
represent the Mott phase, as we have already stated in 
the introduction. 

The fact that there exists no Gaussian state corre- 
sponding to the Mott phase except the state of the an- 
tiferromagentic phase also leads to an incorrect transi- 
tion temperature from the antiferromagnetic to the nor- 
mal phase. For the exact solution of the Hubbard model 
the antiferomagnetic correlations are predicted to be de- 
stroyed at a temperature Tc ~ t^ /u which is the energy 
scale of the superexchange energy. However, in our ap- 
proach the antiferromagnetic correlations appear as soon 
as we find one particle per site, leading to a transition 
temperature Tc ^ u (cf. Fig 14). 



V. THE HUBBARD MODEL IN AN EXTERNAL 
TRAP 

While the translationally invariant Hubbard model al- 
lows for an elegant analytic solution within generalized 
Hartree Fock theory, experimental setups break this sym- 
metry. Since the particles have to be spatially confined, 
an external trapping potential is needed. In many setups, 
this trapping is realized via a harmonic confinement. The 
Hamiltonian of the system is then altered by a position- 
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FIG. 15: Density at the center of the trap for the ground 
state as a function of the chemical potential fi and the in- 
teraction strength u for three different trapping potentials 
Vt = 0.05, 0.1, 0.25 (b) - (d), as well as for the translationally 
invariant case (a). 



dependent chemical potential, i.e. 



{x,y)eA,a 



^ ^ rixt^xi + ^ l^xrix^a, (30) 



ieA 



where for x = x{h^ v) with coordinates (/i, v) we have 



l^x = I^^Vt 



iv^_,y /iv^_;^- 



• (31) 



In the following we investigate how the external confine- 
ment alters the physical properties of the attractive Hub- 
bard model (with open boundary conditions) at zero and 
finite temperature. To gain first insights into the be- 



havior of the system, we have depicted in Fig. [15] the 
density at the center of the trap as a function of the in- 
teraction strength and the chemical potential for three 
different trapping potentials (Fig. IlSlb-d), as well as for 
the translationally invariant case (Fig.flSla). We see that 
the filling in the center of the trap is only altered slightly 
be the external potential. In the following, we call a sys- 
tem half-filled if there is one particle at the center of the 
trap. We find that the condition for half- filling does not 
depend much on the trapping potential, and we can use 
the results from Fig. [15] to adjust the chemical potential 
in order to obtain half-filling in the following. 



A. Attractive Hubbard model in a trap 

We start with a closer investigation of the effects of an 
external trapping potential on the physics of the attrac- 
tive Hubbard model. In Fig. [16] we present the density 
distribution for ix = — 5 at half- filling for four different 



trapping potentials Vt = 0.01, 0.05, 0.1, 0.25. We see that 
for Vt = 0.01 the density distribution is close to that 
of a free system with open boundary conditions, so that 
boundary effects are expected in this case. Hence, we 
will consider only Vt = 0.05,0.1,0.25 from now on, and 
compare the results to the translationally invariant model 
with periodic boundary conditions. 
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FIG. 16: Density profiles for the ground state ( T = 0) for 
the four different trapping potentials Vt = 0.01,0.05,0.1,0.25 
at half-filling (i.e. one particle at the center of the trap) and 
u = —5. 

Now we investigate how the phase diagram of the pair- 
ing changes under the influence of an external confine- 



ment. In Fig 17 we have depicted how the pairing per 
particle as a function of the interaction strength u and 
the chemical potential fi changes under the influence of 
an external trap. We have added lines for 50, 100 and 
150 particles as a guide to the eye. 

To gain further insight into the influence of the trap. 



we have depicted in Fig. [18] the pairing per particle for a 
half-filled trap, i.e. a setting where we find one particle 
at the center of the trap, and compare with the trans- 
lationally invariant case. Note that for a translationally 
invariant system an additional Gauge symmetry emerges 
at half- filling [11], and the underlying Gauge transforma- 
tion does not leave the pairing invariant. However, this 
Gauge symmetry is broken for an inhomogeneous sys- 
tem, and we can only perform a qualitative comparison 
of the translationally invariant and trapped system. In 
this respect, we observe that the pairing decreases with 
decreasing interaction strength, both in the translation- 
ally invariant and in the trapped system, as expected. 
Further, we see that a very strong confinement can lead 
to a slight decrease in the pairing per particle. This ef- 
fect can be understood in the following way: Since we 
fix the number of particles at the center of the trap to 
one, and the external trap acts like a position-dependent 
chemical potential, a strong trap will significantly reduce 
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FIG. 17: Pairing of the ground state ( T = 0) of the at- 
tractive Hubbard model with an external confinement Vt = 
0.05, 0.1, 0.25, as well as for the translationally invariant case 
{Vt = 0), for an interaction range of u ^ [—10,0] and for a 
chemical potential /i G [0, 5]. 




FIG. 18: Pairing per particle for the ground state (T — 0) for 
the attractive Hubbard model with an external confinement 
Vt — 0.05,0.1,0.25 at half-filling, as well as for the transla- 
tionally invariant case [Vt — 0). (Note that in that case the 
absolute value of the pairing is not uniquely determined due 
to an additional Gauge freedom - see text). If the external 
trapping potential becomes too strong, the paired phase is 
suppressed. 



the number of particles at the boundary of the lattice, 
and the paired phase is suppressed. 

We continue our investigations on the influence 
of the external trapping potential by analyzing the 
temperature-dependence of the pairing in an external 
trap. In Fig. 19 we consider interactions i^ = — 10 and 
—5 and a temperature regime of /3 G [0,2]. First, we 



observe that the external trapping alters only slightly 
the temperature-dependence of the pairing. Further, the 
value of /3c, where a phase transition from a paired to a 
normal state occurs is not altered by the external con- 
finement, and is identical to the translationally invariant 
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FIG. 19: Influence of an external trap on the pairing as a func- 
tion of the temperature for i^ = — 10 and i^ = — 5 and trapping 
potentials Vt = 0.05,0.1,0.25, as well as for the translation- 
ally invariant case (Vt = 0) at half-filling. We find that the 
external confinement has no influence on the temperature /3c 
where the transition from the paired to the unpaired phase 
occurs. Note that if the temperature is high enough, the gHF 
Gibbs state becomes unique for the translationally invariant 
case (see text). 



case. At this point we would also like to mention that 
the Gauge freedom of the ground state is lifted with in- 
creasing temperature, and a unique gHF ground state is 
expected if the temperature becomes high enough [11] . 

we can now con- 



With the results depicted in Fig. 19 



elude that for high enough temperature the pairing as 
well as the transition from a paired to an unpaired phase 
in a trap are also approximated well quantitatively by a 
translationally invariant system. 

Finally, we address another question related to super- 
conductivity, namely the symmetry of the pair wave func- 
tion, i.e. the question if we have an s- or p-wave super- 
conductor. Every pure Gaussian state has a wave func- 
tion of the form |^) = exp E:r,^^^,^'(^' ^) 4,^4,(7' l^)- 
Since we consider a translationally invariant system, we 



have (j)a,a'{x^y) = (j)cr^a'{^)^ where f = x — y. In Fig. 20 
we have depicted the pair wave function in momentum 
space, |^t,i(^i'^2)P for a translationally invariant and 
trapped system foiu = —6 and half- filling. We find 
for all trapping potentials a spherically symmetric wave 
function in agreement with s-wave superconductivity. 

In summary, we find that a weak external confinement 
alters only slightly the physics of ground and thermal 
state of the attractive Hubbard model, so that experi- 
ments that require only a weak confinement will be well 
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FIG. 20: Pair wave function in momentum space, |(/)(A:i, /c2)|^ 
for the translationally invariant (Vt = 0) and trapped {Vt — 
0.05,0,1,0.25) system for n = — 6 at half-filling. The wave- 
function is always spherically symmetric, in agreement with 
s-wave superconductivity. 



suited to gain further understanding of the 2d Hubbard 
model. 



VI. DYNAMICAL PROPERTIES OF THE 
HUBBARD MODEL 



FIG. 21: BEC-BCS- crossover for the translationally invariant 
Hubbard model at half-filling. Starting from the ground state 
for u — —7 we perform a linear ramp to the final value u — 7 
with ramping times T/ = 1, 10, 100 and compare the energy 
during the evolution with the ground state energy for various 
u (black triangles) . The evolution is adiabatic for all negative 
u when T/ = 100, while adiabaticity is no longer given for u 
positive. This problem cannot be resolved by longer ramping 
times (a). Further, starting from u > ^ and ramping to i/ < 
leads to the same observation: As soon as i^ = is reached 
the evolution of the system does not follow the ground state 
energy any more (b). 



Current experimental setups with ultra cold gases in 
optical lattices have opened the door to an investigation 
of the dynamical behavior of fermionic lattice systems. 
To this end, it has become possible to study processes 
like the BEC-BCS-crossover or dynamical quenches in 
the laboratory. Unfortunately, existing numerical tech- 
niques, like QMC, are not capable of describing the 
physics of time-evolution, so that numerical benchmark 
results for these processes are still missing so far. 

In this Section we apply the machinery of time- 
evolution within gHFT derived in Sec. |III A| to dynamical 
processes as they can be realized by current experimen- 
tal setups. In the following we consider both, the case of 
a translationally invariant system with periodic bound- 
ary conditions as well as the case of an external trapping 
potential. 



to the final value u = 7. We have used ramping times 
Tf = 1,10,100,200,400 and have compared the energy 
of the evolving system, E{u^Tf)^ with its ground state 
energy, Eq{u)^ for a given u. 

As a result, we observe first that E{u,Tf) has the same 
form for Tj = 100, 200, 400, so that we have only depicted 
the evolutions for Tj = 1, 10, 100 in Fig.[2l] We find that 
for Tf = 100 the evolution is adiabatic for all negative 
u. However, from u = onwards, the energy starts to 
deviate more and more from the true ground state energy. 
An increase in ramping time does not help to resolve this 
problem (Fig.[21p.). In Fig. 21 3 we have depicted a ramp 
from u = 2 to u = —2, and also find that adiabaticity is 
destroyed when we cross the value u = 0. 



Dynamical process of the translationally 
invariant Hubbard model 



In this Section we want to learn more about dynam- 
ical processes of the translationally invariant Hubbard 
model with periodic boundary conditions. To this end 
we have investigated a linear ramp of the interaction 
strength u for various ramping times Tj, addressing 
the question if we can achieve an adiabatic evolution of 
the system. As an example, we have started from the 
ground state for i^ = — 7 and have increased u linearly 



We can gain an understanding of this phenomena by 
looking at Fig. [18] We see that for i^ = the ground 
state undergoes a phase transition from a paired to an 
unpaired phase. However, as we have depicted in Fig. [22} 
the dynamical evolution of the system does not show this 
transition: When we start from li = — 2, where we are 
in a paired phase, the system remains paired even when 
we ramp to the repulsive side (solid lines). On the other 
hand, the unpaired system from which we start for i^ = 2 
remains unpaired even when we ramp to positive values 
of u (dashed lines). 
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FIG. 22: Evolution of the pairing per particle during a linear 
ramp. The solid lines indicate the evolution of the pairing 
under a linear ramp from u — —2 to u — 2. We see that 
the pairing survives even on the repulsive side. Conversely, 
when we start in the unpaired regime ai u = 2 and ramp to 
negative u the pairing remains zero (dashed lines). 



Dynamical process for the Hubbard model in a 
trap 



Finally, we investigate dynamical properties of the 
Hubbard model in a trap. Here, we want to investigate 
how the attractive Hubbard model reacts to a change of 
the trapping potential Vf. To this end, we start with the 
ground state of the system in a shallow trap of strength 
Vt = 0.1 and then squeeze the trap, by changing the trap- 
ping potential linearly in time, to a final trapping poten- 
tial of Vt = 0.25. As an example, we have considered 
an interaction strength of u = —6 and we have set the 
chemical potential to half- filling (here: /i = 3) (I). When 
we have reached the final value Vt = 0.25 we reverse the 
process and now decrease the trapping potential linearly 
to its initial value of Vt = 0.1 (H). 

As the first question, we study if our algorithm allows 
for an adiabatic evolution, i.e. the energy of the system 
for process (I) and (H) is the same. To this end we have 
investigated ramping times of Tj = 100 and Tj = 200, 



and find a reversible process in both cases (see Fig. 23). 



Next, we are interested in the evolution of the pair- 
ing under the change of the trapping potential. The 
result is depicted in Fig. l24| We find that the process 



gets reversible when we go to longer ramping times, and 
that the evolution of the pairing gets smoother. Further, 
we can understand the decrease (increase) of the pairing 
strength with an increasing (decreasing) trapping poten- 
tial by looking at Fig. [6J The change of the trapping 
potential can be understood as a change of the chemical 
potential ja for fixed value of u. 
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FIG. 23: Evolution of the energy for the half- filled Hubbard 
model with i^ = — 6, /i = 3 (half- filling) under a change of 
the external trapping potential for ramping times T/ = 100 
(left) and T/ = 200 (right). Starting from the ground state 
for 14 = 0.1 we change the trapping potential Vt linearly to its 
final value Vt = 0.25 (blue line). Then we reverse the process, 
decreasing Vt again to Vt = 0.1 (purple line). We see that the 
process is reversible. 
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FIG. 24: Evolution of the pairing for the half- filled Hubbard 
model with u = —6, /i = 3 (i.e. half-filling) under a change of 
the external trapping potential for ramping times Tf = 100 
(left) and Tf = 200 (right). Starting from the ground state 
for 14 = 0.1 we change the trapping potential Vt linearly to its 
final value Vt = 0.25 (blue line). Then we reverse the process, 
decreasing Vt again to 14 = 0.1 (purple line). We see that the 
process gets reversible when we go to longer ramping times, 
and that the evolution of the pairing gets smoother. 



Vn. CONCLUSION 

In this work we have formulated the gHFT for 
fermionic lattice systems with a two-body interaction in 
terms of the covariance matrix, and have derived numer- 
ical methods to find the time-evolution, the ground and 
thermal states in this framework. We have benchmarked 
our methods with the two-dimensional translationally in- 
variant Hubbard model, and have found that our numer- 
ical methods work very well to solve the general Hartree- 
Fock theory. Our algorithms run fast, stable and can be 
applied to inhomogeneous systems, as well as to large 
system sizes, since all algorithms are formulated in terms 
of the CM, leading to a quadratic scaling in the system 
size. However, we also found certain limitations of our 
approach that are due to the fact that gHFT is intrinsi- 
cally not capable of describing certain strongly correlated 
phases, like the Mott phase at finite temperature. 

In this respect, we believe that the numerical meth- 
ods developed in this work represent a powerful tool to 
simulate fermionic systems in optical lattices, at least in 
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the weakly interacting regime. Further, we have demon- 
strated that our algorithms can be easily applied to in- 
homogeneous systems, as they are realized, e.g., by the 
Hubbard model with an external trapping potential. 

Concerning the question if properties of the system in 
the thermodynamic limit can be obtained via finite-size 
scaling the results presented in Fig. [7| and table [l] suggest 
that already a moderate system size of about 100 sites 
is enough to get a good estimate of the properties of the 
attractive Hubbard model in the thermodynamic limit. 
In addition, our numerics done for a lattice size up to 
16 X 16 show that energies can be obtained with only 
a slight increase in the running time, while correlation 
functions need some more effort, but can still be obtained 
in a reasonable amount of time 

Thus, in a future work, we aim at applying our ap- 
proach to some particular experimental setup, which will 
require an extension of our current algorithms to much 
larger system sizes, and constitute a further benchmark 
of our methods. 
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Appendix A: Minimization of the energy 

In order to obtain the generalized Hartree-Fock ground 
state we have to solve the minimization problem 



min E{p) = min tT[Hp] 

p Cjaussian p Gaussian 



3 'i'J^k,l 



The condition zF < 1 is fulfilled iff 1 -zF = {A^iB){A^ 
iB)\ where A and B are real matrices [42^, from which 
we derive 



F = AB^ -BA^, 
1 = AA^ ^BB^. 



(Al) 
(A2) 



Then, we can reformulate the problem as an optimization 
problem using Lagrange multipliers. The Lagrangian is 
given by 

3 ^ Uijki{AB'^ - BA^)ij{AB^ - BA^)ki 

ijkl 

+ 2_^ ^kii^ki^ii + BkiBii - Ski), (A3) 

i,k,l 



where A-^ = A is the matrix of the Lagrange multipliers. 
The necessary conditions for a local minimum are given 

by 



OCe 
dAap 
dCE 



dB, 



aP 



2{hB + XA)ap = 0, (A4) 

2{-hA + \B)afs = 0. (A5) 



From this we can obtain equations for the CM if we use 
Eqs. (TaTI) and (|A2|: 



hT + \l = 0, 
h + \T = 0. 



(A6) 
(A7) 



Since F^ = —1 and A-^ = A the condition [h.,T] = 
follows immediately. 



Appendix B: Imaginary time evolution of the CM 

In this Section we give some details that lead to 
Eq. (19). To this end, we will show that either an evo- 



lution with the quadratic, state-dependent Hamiltonian 
Hq or the evolution with the interaction Hamiltonian 
H{T^ U) followed by an application of Wick's theorem 
lead to the desired equation. 



1. Imaginary time evolution with the quadratic, 
state-dependent Hamiltonian 

In the following, we will show how the CM evolves 
under an evolution in imaginary time with the quadratic, 
state-dependent Hamiltonian Hq = i ^- • hijCiCj^ i.e., we 
have to calculate (c.f. Eq. ([iS])) 



ki 



-it! [{Hq, CkCi}] + 2tT[HQp]Tkl. 



First, note that tT[Hp] = — tr[F/i]. Next, we want to 
calculate — ztr [{i^fg, c/cq}] using Wick's theorem: 

- itT[p{t){HQ,CkCi}] = ^hij{ciCjCkCi -^CkCiCiCj) 
= 2 2_^ hij{—rijTki + 2TikTji) — Ahki 

ijy^k,l 

= 2tT[hr]rki-^rhT^h)ki, (bi) 

and we obtain f = -4{ThT + h). 



2. Evolution with the interaction Hamiltonian 

In the following, we will show how the CM evolves 
under an evolution in imaginary time with the in- 
teraction Hamiltonian H{T,U) = ^^ijTijCiCj + 
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Hijki^ijkiCiCjCkCi, i.e. we have to evaluate (c.f. 
Eq. PI) 



tki = -itr [{ff (T, [/), CkCi}] + 2tr[iJ(r, U)p\Tki. 

We split the Hamiltonian H(T^ U) in two terms, Hq = 
'^Y.ki'^kiCkCi and Hi = ^j^i^^UkimnCkCiCmCn, and per- 
form the calculation for both terms independently. First, 
the calculation for Hg is t he same as for Hq, so that we 
can deduce from Eq. (iBll) 



-ztr[p(t){i7„CfcQ}] =2tr[T^]^fc,-4(^T^ + T)fc^ (B2) 

Next, we calculate the contribution of the interaction 
term. First, 

-itT[{Hi,CkCi}p{t)] = 

-2i ^ UijmniCiCjCmCnCkCi) - 24tTB[UT]kl. 

i,j,'m,n^k,l 

Next, using Wick's theorem, we find 

^2- ^ ^ ^ijiTinxCiCjCfYiCfiCj^Ci) = 

i,j,m,ny^k,l 

/ ^ ^ijmn Vt)i ijL rnn^ kl ^4i ij^L jm^ nl) • 

i,j,'m,n^k,l 

We treat the two terms independently, starting with 

Z^i,j,m,ny^k,l ^ijmn^ ij^ mn- 

/ ^ ^ijmn^ ij^ ran ^ 

tr[rtrB[C/r]] - 4(trB[C/r]r)fcfc - A{ivB[UT]T)ii 

- 8 ^ UijkiTikTji - Ail B[UT]kiTki. 

Next, we calculate Y,i,j,m,n^k,l UijmnTikT'jniTnl- 

/ ^ ^ijmn^ i/e-L jm^ nl — 

i,j,m,nyi^k,l 

{rtTB[ur]r)ki - {tTBiUTjr)^^^^ - {tTB[UT]r)urki 

— tTB[UT]klTf^l — 2 ^ ^ Uklin^ik^nl^kl 



Thus, we obtain 

-itT[p{Hi,CkCi}] = 

- 24tT[UT]ki - 24(rtr[/7r]r)fez + 6tr[rtr[/7r]]rfcz. 

Hence, -itT[p{H,CkCj\] = -2tT[H p]T m - ^ThT ^ h) , 
which leads to Eq. ( [T9| . 



Appendix C: Minimization of the free energy 



We want to minimize 

min F{p)= min {E{p) - r'S{p)} . 

p Gaussian p Gaussian 



with the help of Lagrange multipliers. An expression 
for the energy has already been given in Eq. (14). The 
entropy S — — tr[plnp] is given by S* = Mln2 — 
|tr [(1 + iF) ln(l + zF)] [43 , so that we have to solve the 
following optimization problem: 

min i^(p) = 

p Gaussian 



mini -tr[(T + 3trB[t/F])F] 
■ /3-M M In 2 - -tr [(1 + iV) ln(l + iV)\ 



\ 



As we have already stated in Appendix |AJ the condition 
zF < 1 is fulfilled iff 1 - zF = (A + iB\A ^ iB)\ where 
A and B are real matrices [42 , from which we derive 



F = AB^ -BA^, 
1 = AA^ ^BB^. 



(CI) 
(C2) 



Thus, we can reformulate the minimization of the free 
energy as an optimization problem with Lagrange multi- 
pliers. The Lagrangian is of the form 

£5= F(A,5) + ^AfcK^fci^H + ^fci^/i-4/), (C3) 

where A = A-^ is the matrix of Lagrange multipliers. 
Now, ^^ = §^w¥^ (resp. B), and we calculate 
first ■§^- To this end, we consider the entropy term, 

tr[(l + iF) ln(l + iV)\. We us ln(l - X) = - YX^i \^^ 
(see e.g. [4T) leading to 

tr[(l + iV) ln(l + iY)\ = - ^ ^(-i)'(r'^ + ir'=+')«. 



kl 



Then 
d 



dV 



kl 



-tr[(l+zF) ln(l+zF)] = z [1 + ln(l - iT)]^^ , (C4) 



where we have made use of the formula (1 — X) ^ = 
E^o^^ [44 . Since [ln(l - iT)f = ln(l + zF), we 
arrive at 



dCs 

^^mn 

PCs 



2{hFB + XA) = 0, 
2{-hFA^XB) =0, 



Then, using Eqs. (CI) and (C2) we obtain the following 
necessary conditions for a minimal entropy: 



[hF{r),r] = 0, 
/if(r)(i + r2) = 0. 



(C5) 
(C6) 
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